% Copyright (C) 2021-2023 Benjamin Born, Francesco D'Ascanio, Gernot J. Mueller, Johannes Pfeifer
%
% This is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% It is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
% 
% For a copy of the GNU General Public License,
% see <http://www.gnu.org/licenses/>.


%Compute asymmetric LP to first-stage VAR shocks identified using sign
%restrictions 

clear
close all
set(0,'defaulttextinterpreter','latex')

if ~isfolder('Figures')
    mkdir('.','Figures');
end

max_irf_horizon = 14;
only_90 = 0;

fontsize = 6;
print_dummy = 1;

%% Data construction
load('results/IRFS_symmetric_1983.mat')
load('results/structural_shocks_point_1983.mat')

VAR_IRFS = PF_G_IRFS;
timeline = timeline';

g_shocks        = structural_shocks.PF_G;
log_g           = YYact(:, 2);
log_y           = YYact(:, 3);
log_rx          = YYact(:, 4);
log_tax         = YYact(:, 1);


%% Empirics
start_date = 1983;
end_date = timeline(end);
% cut sample
g_shocks       = g_shocks(timeline>=start_date & timeline<=end_date);
log_g        = log_g(timeline>=start_date & timeline<=end_date);
log_y  = log_y(timeline>=start_date & timeline<=end_date);
log_rx        = log_rx(timeline>=start_date & timeline<=end_date);
log_tax       = log_tax(timeline>=start_date & timeline<=end_date);

n_obs = size(log_g, 1);
t = (1:n_obs)';
t2 = t.^2 / sqrt(n_obs);

dependent_var_names = {'Government spending', 'Output', 'Real Exchange Rate', 'Tax revenues'};
dependent_var_mat = [log_g, log_y, log_rx, log_tax];
xlabel_vec = {'percent','percent','percent','percent'};
dependent_var_mat(any(isnan(dependent_var_mat), 2), :) = NaN;

%% Asymmetric

g_shocks_negative = zeros(size(g_shocks));
g_shocks_negative(g_shocks<0 | isnan(g_shocks)) = g_shocks(g_shocks<0 | isnan(g_shocks));
g_shocks_positive = zeros(size(g_shocks));
g_shocks_positive(g_shocks>0 | isnan(g_shocks)) = g_shocks(g_shocks>0 | isnan(g_shocks));
if any(any(g_shocks_negative>0)) || any(any(g_shocks_positive<0 ))
    error('Sign is wrong')
end

regressor_base = [g_shocks_negative g_shocks_positive t lagmatrix([log_g, log_y, log_rx], [1 2 3 4])];

maxNWLag = floor(4*(n_obs/100)^(2/9)); % for Newey-West
theta_mat = NaN(1+max_irf_horizon, length(dependent_var_names), 2);
se_mat = NaN(1+max_irf_horizon, length(dependent_var_names), 2);

for var_iter=1:length(dependent_var_names)
    resid = [];
    for horizon_iter = 0:max_irf_horizon
        regressor_mat_temp = [regressor_base resid];
        dep_var_temp = lagmatrix(dependent_var_mat(:, var_iter), -horizon_iter);

        [~, se, coeff] = hac(regressor_mat_temp, dep_var_temp,'bandwidth', maxNWLag+1,'display','off'); % Newey-West regression
        resid = dep_var_temp - [ones(n_obs, 1) regressor_mat_temp] * coeff;
        theta_mat(1+horizon_iter, var_iter, 1) = coeff(2);
        se_mat(1+horizon_iter, var_iter, 1)  = se(2);
        theta_mat(1+horizon_iter, var_iter, 2) = coeff(3);
        se_mat(1+horizon_iter, var_iter, 2)  = se(3);
    end
end

%% Plotting

confidence_90  = norminv(0.9 + (1 - 0.9) / 2, 0, 1);
confidence_68  = norminv(0.68 + (1 - 0.68) / 2, 0, 1);
mycolors = [1 1 1; .85 .85 .85];

main_fig = figure;
colormap(mycolors);

for var_iter=1:length(dependent_var_names)
    subplot(4,2,(var_iter-1)*2+1)
    if var_iter ==3
        scaling = 100/(g_shocks\log_g*100);
    else
        scaling = -100/(g_shocks\log_g*100);
    end
    ci1_68 = scaling*(theta_mat(:, var_iter, 1)+confidence_68*se_mat(:, var_iter, 1));
    ci2_68 = scaling*(theta_mat(:, var_iter, 1)-confidence_68*se_mat(:, var_iter, 1));
    topci_68 = max(ci1_68,ci2_68);
    bottomci_68 = min(ci1_68,ci2_68);
    ci1_90 = scaling*(theta_mat(:, var_iter, 1)+confidence_90*se_mat(:, var_iter, 1));
    ci2_90 = scaling*(theta_mat(:, var_iter, 1)-confidence_90*se_mat(:, var_iter, 1));
    topci_90 = max(ci1_90,ci2_90);
    bottomci_90 = min(ci1_90,ci2_90);
    ha1_90 = area(0:max_irf_horizon,[bottomci_90, topci_90-bottomci_90],'FaceColor',[204/255 229/255 1],'EdgeColor','none','ShowBaseLine','off'); %[204/255 229/255 1]
    set(ha1_90(1), 'FaceColor', 'none') % this makes the bottom area invisible
    set(ha1_90, 'LineStyle', '-')
    hold on
    if ~only_90
        ha1_68 = area(0:max_irf_horizon,[bottomci_68, topci_68-bottomci_68],'FaceColor',[153/255 204/255 1],'EdgeColor','none','ShowBaseLine','off');
        set(ha1_68(1), 'FaceColor', 'none') % this makes the bottom area invisible
        set(ha1_68, 'LineStyle', '-')
        hold on
    end
    hold on
    fs=plot(0:max_irf_horizon,scaling*theta_mat(:, var_iter, 1),'b-', 'LineWidth', 2.5);
    hline(0,'k:')
    xlim([0 max_irf_horizon])
    box on;set(gca,'xTick',0:2:max_irf_horizon,'Layer','top','FontSize',fontsize);
    xtickangle(0)
    title(dependent_var_names(var_iter),'FontSize',fontsize)
    set(0,'DefaultAxesTitleFontWeight','normal');
    set(gca,'TickLabelInterpreter','latex')
    ylabel(xlabel_vec{var_iter},'FontSize', fontsize)
    if var_iter==3
        xlabel('quarters','FontSize', fontsize)
    end
    set(gca, 'FontSize', fontsize)

    subplot(4,2,2*var_iter)
    if var_iter ==3
        scaling = -100/(g_shocks\log_g*100);
    else
        scaling = 100/(g_shocks\log_g*100);
    end
    ci1_68 = scaling*(theta_mat(:, var_iter, 2)+confidence_68*se_mat(:, var_iter, 2));
    ci2_68 = scaling*(theta_mat(:, var_iter, 2)-confidence_68*se_mat(:, var_iter, 2));
    topci_68 = max(ci1_68,ci2_68);
    bottomci_68 = min(ci1_68,ci2_68);
    ci1_90 = scaling*(theta_mat(:, var_iter, 2)+confidence_90*se_mat(:, var_iter, 2));
    ci2_90 = scaling*(theta_mat(:, var_iter, 2)-confidence_90*se_mat(:, var_iter, 2));
    topci_90 = max(ci1_90,ci2_90);
    bottomci_90 = min(ci1_90,ci2_90);
    ha1_90 = area(0:max_irf_horizon,[bottomci_90, topci_90-bottomci_90],'FaceColor',[204/255 229/255 1],'EdgeColor','none','ShowBaseLine','off'); %[204/255 229/255 1]
    set(ha1_90(1), 'FaceColor', 'none') % this makes the bottom area invisible
    set(ha1_90, 'LineStyle', '-')
    hold on
    if ~only_90
        ha1_68 = area(0:max_irf_horizon,[bottomci_68, topci_68-bottomci_68],'FaceColor',[153/255 204/255 1],'EdgeColor','none','ShowBaseLine','off');
        set(ha1_68(1), 'FaceColor', 'none') % this makes the bottom area invisible
        set(ha1_68, 'LineStyle', '-')
        hold on
    end
    hold on
    fs=plot(0:max_irf_horizon,scaling*theta_mat(:, var_iter, 2),'b-', 'LineWidth', 2.5);
    hline(0,'k:')
    xlim([0 max_irf_horizon])
    box on;set(gca,'xTick',0:2:max_irf_horizon,'Layer','top','FontSize',fontsize);
    xtickangle(0)
    title(dependent_var_names(var_iter),'FontSize',fontsize)
    set(0,'DefaultAxesTitleFontWeight','normal');
    set(gca,'TickLabelInterpreter','latex')
    ylabel(xlabel_vec{var_iter},'FontSize', fontsize)
    if var_iter==3
        xlabel('quarters','FontSize', fontsize)
    end
    set(gca, 'FontSize', fontsize)
end
saveas(main_fig,'Figures/us_asym_signrestrict');

%% Symmetric
regressor_base = [g_shocks lagmatrix([log_g, log_y, log_rx, log_tax], [1 2 3 4])];

maxNWLag = floor(4*(n_obs/100)^(2/9)); % for Newey-West
theta_mat = NaN(1+max_irf_horizon, length(dependent_var_names));
se_mat = NaN(1+max_irf_horizon, length(dependent_var_names));

for var_iter=1:length(dependent_var_names)
    resid = [];
    for horizon_iter = 0:max_irf_horizon
        regressor_mat_temp = [regressor_base resid];
        dep_var_temp = lagmatrix(dependent_var_mat(:, var_iter), -horizon_iter);
        [~, se, coeff] = hac(regressor_mat_temp, dep_var_temp,'bandwidth', maxNWLag+1,'display','off'); % Newey-West regression
        resid = dep_var_temp - [ones(n_obs, 1) regressor_mat_temp] * coeff;
        theta_mat(1+horizon_iter, var_iter) = coeff(2);
        se_mat(1+horizon_iter, var_iter)  = se(2);
    end
end

%% Plotting

confidence_90  = norminv(0.9 + (1 - 0.9) / 2, 0, 1);
confidence_68  = norminv(0.68 + (1 - 0.68) / 2, 0, 1);
mycolors = [1 1 1; .85 .85 .85];

main_fig = openfig('Figures/us_asym_signrestrict');

colormap(mycolors);

for var_iter=1:length(dependent_var_names)
    subplot(4,2,(var_iter-1)*2+1)
    scaling = -1;
    ci1_68 = scaling * VAR_IRFS(1:max_irf_horizon+1,var_iter,2);
    ci2_68 = scaling * VAR_IRFS(1:max_irf_horizon+1,var_iter,4);
    topci_68 = max(ci1_68,ci2_68);
    bottomci_68 = min(ci1_68,ci2_68);
    ci1_90 = scaling * VAR_IRFS(1:max_irf_horizon+1,var_iter,1);
    ci2_90 = scaling * VAR_IRFS(1:max_irf_horizon+1,var_iter,5);
    topci_90 = max(ci1_90,ci2_90);
    bottomci_90 = min(ci1_90,ci2_90);
    ha1_90 = area(0:max_irf_horizon,[bottomci_90, topci_90-bottomci_90],'FaceColor',[1 204/255 204/255],'EdgeColor','none','ShowBaseLine','off'); %[204/255 229/255 1]
    set(ha1_90(1), 'FaceColor', 'none') % this makes the bottom area invisible
    set(ha1_90, 'LineStyle', '-')
    ha1_90(2).FaceAlpha = 0.5;
    hold on
    if ~only_90
        ha1_68 = area(0:max_irf_horizon,[bottomci_68, topci_68-bottomci_68],'FaceColor',[1 0.6 0.6],'EdgeColor','none','ShowBaseLine','off');
        set(ha1_68(1), 'FaceColor', 'none') % this makes the bottom area invisible
        set(ha1_68, 'LineStyle', '-')
        ha1_68(2).FaceAlpha = 0.5;
        hold on
    end
    fs=plot(0:max_irf_horizon,scaling * VAR_IRFS(1:max_irf_horizon+1,var_iter,3),'r--', 'LineWidth', 2.5);
    hline(0,'k:')
    xlim([0 max_irf_horizon])
    box on;set(gca,'xTick',0:2:max_irf_horizon,'Layer','top','FontSize',fontsize);
    xtickangle(0)
    title(dependent_var_names(var_iter),'FontSize',fontsize)
    set(0,'DefaultAxesTitleFontWeight','normal');
    set(gca,'TickLabelInterpreter','latex')
    ylabel(xlabel_vec{var_iter},'FontSize', fontsize)
    if var_iter==3
        xlabel('quarters','FontSize', fontsize)
    end
    set(gca, 'FontSize', fontsize)

    subplot(4,2,2*var_iter)
    scaling = 1;
    ci1_68 = scaling * VAR_IRFS(1:max_irf_horizon+1,var_iter,2);
    ci2_68 = scaling * VAR_IRFS(1:max_irf_horizon+1,var_iter,4);
    topci_68 = max(ci1_68,ci2_68);
    bottomci_68 = min(ci1_68,ci2_68);
    ci1_90 = scaling * VAR_IRFS(1:max_irf_horizon+1,var_iter,1);
    ci2_90 = scaling * VAR_IRFS(1:max_irf_horizon+1,var_iter,5);
    topci_90 = max(ci1_90,ci2_90);
    bottomci_90 = min(ci1_90,ci2_90);
    ha1_90 = area(0:max_irf_horizon,[bottomci_90, topci_90-bottomci_90],'FaceColor',[1 204/255 204/255],'EdgeColor','none','ShowBaseLine','off'); %[204/255 229/255 1]
    set(ha1_90(1), 'FaceColor', 'none') % this makes the bottom area invisible
    set(ha1_90, 'LineStyle', '-')
    ha1_90(2).FaceAlpha = 0.5;
    hold on
    if ~only_90
        ha1_68 = area(0:max_irf_horizon,[bottomci_68, topci_68-bottomci_68],'FaceColor',[1 0.6 0.6],'EdgeColor','none','ShowBaseLine','off');
        set(ha1_68(1), 'FaceColor', 'none') % this makes the bottom area invisible
        set(ha1_68, 'LineStyle', '-')
        ha1_68(2).FaceAlpha = 0.5;
        hold on
    end
    fs=plot(0:max_irf_horizon,scaling * VAR_IRFS(1:max_irf_horizon+1,var_iter,3),'r--', 'LineWidth', 2.5);
    hline(0,'k:')
    xlim([0 max_irf_horizon])
    box on;set(gca,'xTick',0:2:max_irf_horizon,'Layer','top','FontSize',fontsize);
    xtickangle(0)
    title(dependent_var_names(var_iter),'FontSize',fontsize)
    set(0,'DefaultAxesTitleFontWeight','normal');
    set(gca,'TickLabelInterpreter','latex')
    ylabel(xlabel_vec{var_iter},'FontSize', fontsize)
    if var_iter==3
        xlabel('quarters','FontSize', fontsize)
    end
    set(gca, 'FontSize', fontsize)
end
saveas(main_fig,'Figures/us_linear_v_asym_signrestrict');
print('Figures/us_linear_v_asym_signrestrict','-dpdf')
